library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.4     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(rgdal)
## Loading required package: sp
## Please note that rgdal will be retired by the end of 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## 
## rgdal: version: 1.5-25, (SVN revision 1143)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.4-5
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
## Overwritten PROJ_LIB was /Library/Frameworks/R.framework/Versions/4.1/Resources/library/rgdal/proj
library(raster)
## 
## Attaching package: 'raster'
## The following object is masked from 'package:plotly':
## 
##     select
## The following object is masked from 'package:dplyr':
## 
##     select
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.2.1, PROJ 7.2.1
library(spatialEco)
## 
## Attaching package: 'spatialEco'
## The following object is masked from 'package:raster':
## 
##     shift
## The following object is masked from 'package:dplyr':
## 
##     combine
library(viridis)
## Loading required package: viridisLite
library(rgdal)
library(spatialreg)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## 
## Attaching package: 'spData'
## The following object is masked from 'package:spatialEco':
## 
##     elev
## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(spgwr)
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
library(raster)

knitr::opts_chunk$set(
  fig.width = 6,
  fig.asp = .6,
  out.width = "90%"
)
theme_set(theme_minimal() + theme(legend.position = "bottom"))
options(
  ggplot2.continuous.colour = "viridis",

  ggplot2.continuous.fill = "viridis"
)
scale_colour_discrete = scale_colour_viridis_d
scale_fill_discrete = scale_fill_viridis_d

#Map cumulative mortality

covid_death_sdf = st_read("deathrate_map.shp") 
## Reading layer `deathrate_map' from data source 
##   `/Users/irene/Desktop/DATA SCIENCE/p8105_finalproject/deathrate_map.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 178 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  GCS_WGS84_DD
fig <- ggplotly(
  ggplot(covid_death_sdf) +
    geom_sf(aes(fill = cum_death)) +
    theme(legend.text = element_text(size = 6), 
        legend.position = 'right', axis.text.x = element_text(angle = 45)) +
    labs(
    title = "Average Death Rate per 100,000 people \n from March - September 2020 by ZCTA",
    x = "Latitude",
    y = "Longitude", fill = "Average Death Rate\n (deaths/100,000)")) 
  
fig

#Map cumulative incidence

covid_death_sdf1 = st_read("deathrate_map1.shp") 
## Reading layer `deathrate_map1' from data source 
##   `/Users/irene/Desktop/DATA SCIENCE/p8105_finalproject/deathrate_map1.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 178 features and 22 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -74.25559 ymin: 40.49612 xmax: -73.70001 ymax: 40.91553
## Geodetic CRS:  WGS 84
fig <- ggplotly(
  ggplot(covid_death_sdf1) +
    geom_sf(aes(fill = incidence)) +
    theme(legend.text = element_text(size = 6), 
        legend.position = 'right', axis.text.x = element_text(angle = 45)) +
    labs(
    title = "Covid-19 cumulative Incidence \n March - September 2020 by ZCTA",
    x = "Latitude",
    y = "Longitude", fill = "Covid-19 cumulative Incidence \n (cases/100,000)")) 
  
fig
fig <- ggplotly(
  ggplot(covid_death_sdf1) +
    geom_sf(aes(fill = pro_offers)) +
    theme(legend.text = element_text(size = 6), 
        legend.position = 'right', axis.text.x = element_text(angle = 45)) +
    labs(
    title = "Proportion of SHSAT offers \n in 2020 by ZCTA",
    x = "Latitude",
    y = "Longitude", fill = "proportion of SHSAT offers")) 
  
fig